home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 August / August CD.bin / Shareware / Education / numericalmethods Folder / chap_7 / a7_5.m < prev    next >
Encoding:
Text File  |  1994-06-05  |  2.3 KB  |  96 lines  |  [MATS/MATL]
  1. echo off;
  2. % NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1994
  3. % To accompany the text:
  4. % NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
  5. % Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
  6. % This free software is complements of the author.
  7.  
  8. % Algorithm 7.5 (Adaptive Quadrature Using Simpson's Rule).
  9. % Section    7.4,    Adaptive Quadrature, Page 389
  10. % This program is to be used for pedagogical purposes.
  11. echo on; clc; format long; hold off; clear
  12.  
  13. % This program implements adaptive quadrature.
  14.  
  15. %    Define and store the function f(x)    in the M-file  f.m 
  16. % function y = f(x)
  17. % y = 13.*(x - x.^2).*exp(-3.*x./2);
  18.  
  19. delete f.m
  20. diary f.m; disp('function y = f(x)');...
  21.            disp('y = 13.*(x - x.^2).*exp(-3.*x./2);');...
  22. diary off;
  23.  
  24. f(0); % Remark. f.m adapt.m srule.m are used for Algorithm 7.5
  25.  
  26. pause    % Press any key to see the graph y = f(x).
  27.  
  28. clc; clg;
  29. a = 0;
  30. b = 4;
  31. h = (b-a)/200;
  32. X = a:h:b;
  33. Y = f(X);
  34. c = -1.5;
  35. d = 2;
  36. axis([a b c d]);...
  37. plot(X,Y,'g');...
  38. hold on;...
  39. plot([a b],[0 0],'b',[0 0],[c d],'b');...
  40. xlabel('x');...
  41. ylabel('y');...
  42. title('The curve y = f(x) over [a,b].');...
  43. grid;...
  44. axis;...
  45. hold off;...
  46. shg; pause    % Press any key to continue.
  47.  
  48. clc;
  49.  
  50. % Adaptive quadrature is applied over the interval [a,b].
  51.  
  52. %    Place the endpoints of [a,b] in  a  and  b.
  53.  
  54. % Place the tolerance in  toler.
  55.  
  56. a  = 0;
  57. b  = 4;
  58. toler = 0.00001;
  59.  
  60. [SRmat,quad,err] = adapt('f',a,b,toler);
  61.  
  62. pause    % Press any key to continue.
  63.  
  64. clc;
  65. Mx1 = 'The adaptive quadrature approximation is: ';
  66. Mx2 = '   quadrature value   +- error bound';
  67. clc,echo off,diary output,...
  68. disp(''),disp(Mx1),disp(''),...
  69. disp(Mx2),disp([quad err]),diary off,echo on
  70.  
  71. pause % Press any key to continue.
  72.  
  73. clg; 
  74. X0 = [SRmat(:,1)',b];
  75. Y0 = f(X0);
  76. Z0 = zeros(1,length(X0));
  77. axis([a b c d]);...
  78. plot(X,Y,'g',X0,Y0,'or',X0,Z0,'+r');...
  79. hold on;...
  80. plot([a b],[0 0],'b',[0 0],[c d],'b');...
  81. xlabel('x');...
  82. ylabel('y');...
  83. title('Adaptive quadrature.');...
  84. grid;...
  85. axis;...
  86. hold off;...
  87. shg; pause    % Press any key to continue.
  88.  
  89. Mx3 = 'Quadrature over each subinterval:';
  90. Mx4 = '   a(k)               b(k)               S(a(k),b(k))';
  91. clc,echo off,diary output,disp(''),disp(Mx3),disp(''),...
  92. disp(Mx4),disp(SRmat(:,[1 2 4])),...
  93. disp(['   ∑ S(a(k),b(k))     ∑ error(k)']),...
  94. disp(Mx2),disp([quad err]),diary off,echo on
  95.  
  96.